---
title: "R Notebook"
output: html_notebook
---

# Preamble + Load Data

```{r}

# install.packages("BiocManager")
# BiocManager::install("multtest")

rm(list = ls())

options(scipen = 99)

require(multtest)
require(rio)
require(tidyverse)

#path = "Datasets"

# Load unadjusted P-Vals and divide them by family of correction 
df.p = import("Datasets/pvals_unadjusted.csv") %>% as_tibble()
colnames(df.p) = colnames(df.p) %>% str_remove('[0-9]')

df.family = df.p %>% 
  filter( var_name %in%  c("earnings", "wellbeing", "cogindex", "pref") ) %>% 
  mutate(order = seq(1,4)) %>% 
  relocate(order)

df.work = df.p %>% 
  filter( var_name %in%  c("prod", "typing", "output"))

df.wb = df.p %>% 
  filter( var_name %in%  c("physical", "psych"))

df.cog = df.p %>% 
  filter( var_name %in%  c("cogfunction", "attention"))

df.pref = df.p %>% 
  filter( var_name %in%  c("timepref", "social", "risk"))

# Load step-down FWER (Anderson 08) adjusted p-vals

df_step_down = import(file.path(path, "pval_adj_table_pool.csv")) %>% as_tibble()

df_step_down_family = df_step_down[1:12,]

# Parameters
# corrections
corrections = c("BH", "BY", "hochberg", "holm", "TSBH")
alpha = 0.05 #(for rejection threshold when using the BKY06 method)

```


Function Adjustment
```{r}

p.adjust.vec = function(p.unadjust, corrections, alpha){
  p.adj.mat = matrix(NA, 
                     ncol = length(corrections), 
                     nrow = length(p.unadjust))
    colnames(p.adj.mat) = corrections
    
  for(i in seq(1, length(corrections))){
    corr = corrections[i]
    
    if(corr == "TSBH"){
      ls.p.adj = multtest::mt.rawp2adjp(p.unadjust,
                                 proc = corr,
                                 alpha = 0.05)  
      p.adj.mat[,i] = ls.p.adj$adjp[,2]
    } else{
      ls.p.adj = p.adjust(p.unadjust, method = corr )
      p.adj.mat[,i] = ls.p.adj
    }
      
    
  }
    
  return(p.adj.mat)
} 


```

# Main table 

```{r}

# Comparison: NS vs. Control

df.family.ns.cont = df.family %>% 
  select(order, var_name, p_val_ns) %>% 
  mutate(comparison = "NS vs. Control",
         order_2 = 1) %>% 
  rename(p_val_unadjust = p_val_ns) %>% 
  arrange(p_val_unadjust) 

p.unadjust = df.family.ns.cont$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.family.ns.cont = df.family.ns.cont %>% 
  bind_cols(p.adj.mat)

# Comparison: Nap vs. Control

df.family.nap.cont = df.family %>%
  select(order, var_name, p_val_nap) %>%
  mutate(comparison = "Nap vs. Control",
         order_2 = 2) %>% 
  rename(p_val_unadjust = p_val_nap) %>% 
  arrange(p_val_unadjust)

p.unadjust = df.family.nap.cont$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.family.nap.cont = df.family.nap.cont %>% 
  bind_cols(p.adj.mat)

# Comparison: Nap vs. NS

df.family.nap.ns = df.family %>% 
  select(order, var_name, p_val_d_nap_ns) %>% 
  mutate(comparison = "Nap vs. NS",
         order_2 = 3) %>% 
  rename(p_val_unadjust = p_val_d_nap_ns) %>% 
  arrange(p_val_unadjust)

p.unadjust = df.family.nap.ns$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.family.nap.ns = df.family.nap.ns %>% 
  bind_cols(p.adj.mat)

df.family.all = df.family.ns.cont %>% 
  bind_rows(df.family.nap.cont) %>% 
  bind_rows(df.family.nap.ns) %>% 
  relocate(var_name, comparison) %>% 
  rename(BKY06 = TSBH) %>% 
  arrange(order_2, order) %>% 
  mutate(ri_step_down = df_step_down_family$p_val_adj) %>% 
  select(!contains("order")) %>% 
  relocate(var_name, comparison, p_val_unadjust, ri_step_down, hochberg, holm, BH, BY, BKY06)

export(df.family.all, "Datasets/pval_adj_family.xlsx")

```

# Work Family

```{r}

# Comparison: NS vs. Control

df.work.ns.cont = df.work %>% 
  select(var_name, p_val_ns) %>% 
  mutate(comparison = "NS vs. Control") %>% 
  rename(p_val_unadjust = p_val_ns) %>% 
  arrange(p_val_unadjust)

p.unadjust = df.work.ns.cont$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.work.ns.cont = df.work.ns.cont %>% 
  bind_cols(p.adj.mat)

# Comparison: Nap vs. Control

df.work.nap.cont = df.work %>% 
  select(var_name, p_val_nap) %>% 
  mutate(comparison = "Nap vs. Control") %>% 
  rename(p_val_unadjust = p_val_nap) %>% 
  arrange(p_val_unadjust)

p.unadjust = df.work.nap.cont$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.work.nap.cont = df.work.nap.cont %>% 
  bind_cols(p.adj.mat)

# Comparison: Nap vs. NS

df.work.nap.ns = df.work %>% 
  select(var_name, p_val_d_nap_ns) %>% 
  mutate(comparison = "Nap vs. NS") %>% 
  rename(p_val_unadjust = p_val_d_nap_ns)  %>% 
  arrange(p_val_unadjust)

p.unadjust = df.work.nap.ns$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.work.nap.ns = df.work.nap.ns %>% 
  bind_cols(p.adj.mat)


df.work.all = df.work.ns.cont %>% 
  bind_rows(df.work.nap.cont) %>% 
  bind_rows(df.work.nap.ns) %>% 
  relocate(var_name, comparison)

df.work.all = df.work.all %>% 
  rename(BKY06 = TSBH) %>% 
  arrange(comparison, var_name)

export(df.work.all, "Datasets/pval_adj_work.xlsx")

```

# Well-Being Family

```{r}

# Comparison: NS vs. Control

df.wb.ns.cont = df.wb %>% 
  select(var_name, p_val_ns) %>% 
  mutate(comparison = "NS vs. Control") %>% 
  rename(p_val_unadjust = p_val_ns) %>% 
  arrange(p_val_unadjust)

p.unadjust = df.wb.ns.cont$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.wb.ns.cont = df.wb.ns.cont %>% 
  bind_cols(p.adj.mat)

# Comparison: Nap vs. Control

df.wb.nap.cont = df.wb %>% 
  select(var_name, p_val_nap) %>% 
  mutate(comparison = "Nap vs. Control") %>% 
  rename(p_val_unadjust = p_val_nap) %>% 
  arrange(p_val_unadjust)

p.unadjust = df.wb.nap.cont$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.wb.nap.cont = df.wb.nap.cont %>% 
  bind_cols(p.adj.mat)

# Comparison: Nap vs. NS

df.wb.nap.ns = df.wb %>% 
  select(var_name, p_val_d_nap_ns) %>% 
  mutate(comparison = "Nap vs. NS") %>% 
  rename(p_val_unadjust = p_val_d_nap_ns)  %>% 
  arrange(p_val_unadjust)

p.unadjust = df.wb.nap.ns$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.wb.nap.ns = df.wb.nap.ns %>% 
  bind_cols(p.adj.mat)


df.wb.all = df.wb.ns.cont %>% 
  bind_rows(df.wb.nap.cont) %>% 
  bind_rows(df.wb.nap.ns) %>% 
  relocate(var_name, comparison)

df.wb.all = df.wb.all %>% 
  rename(BKY06 = TSBH) %>% 
  arrange(comparison, var_name)

export(df.wb.all,  "Datasets/pval_adj_wb.xlsx")

```

# Cognition Family

```{r}

# Comparison: NS vs. Control

df.cog.ns.cont = df.cog %>% 
  select(var_name, p_val_ns) %>% 
  mutate(comparison = "NS vs. Control") %>% 
  rename(p_val_unadjust = p_val_ns) %>% 
  arrange(p_val_unadjust)

p.unadjust = df.cog.ns.cont$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.cog.ns.cont = df.cog.ns.cont %>% 
  bind_cols(p.adj.mat)

# Comparison: Nap vs. Control

df.cog.nap.cont = df.cog %>% 
  select(var_name, p_val_nap) %>% 
  mutate(comparison = "Nap vs. Control") %>% 
  rename(p_val_unadjust = p_val_nap) %>% 
  arrange(p_val_unadjust)

p.unadjust = df.cog.nap.cont$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.cog.nap.cont = df.cog.nap.cont %>% 
  bind_cols(p.adj.mat)

# Comparison: Nap vs. NS

df.cog.nap.ns = df.cog %>% 
  select(var_name, p_val_d_nap_ns) %>% 
  mutate(comparison = "Nap vs. NS") %>% 
  rename(p_val_unadjust = p_val_d_nap_ns) %>% 
  arrange(p_val_unadjust)

p.unadjust = df.cog.nap.ns$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.cog.nap.ns = df.cog.nap.ns %>% 
  bind_cols(p.adj.mat)


df.cog.all = df.cog.ns.cont %>% 
  bind_rows(df.cog.nap.cont) %>% 
  bind_rows(df.cog.nap.ns) %>% 
  relocate(var_name, comparison)

df.cog.all = df.cog.all %>% 
  rename(BKY06 = TSBH) %>% 
  arrange(comparison, var_name)

export(df.cog.all,  "Datasets/pval_adj_cog.xlsx")

```

# Preferences Family

```{r}

# Comparison: NS vs. Control

df.pref.ns.cont = df.pref %>% 
  select(var_name, p_val_ns) %>% 
  mutate(comparison = "NS vs. Control") %>% 
  rename(p_val_unadjust = p_val_ns) %>% 
  arrange(p_val_unadjust)

p.unadjust = df.pref.ns.cont$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.pref.ns.cont = df.pref.ns.cont %>% 
  bind_cols(p.adj.mat)

# Comparison: Nap vs. Control

df.pref.nap.cont = df.pref %>% 
  select(var_name, p_val_nap) %>% 
  mutate(comparison = "Nap vs. Control") %>% 
  rename(p_val_unadjust = p_val_nap) %>% 
  arrange(p_val_unadjust)

p.unadjust = df.pref.nap.cont$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.pref.nap.cont = df.pref.nap.cont %>% 
  bind_cols(p.adj.mat)

# Comparison: Nap vs. NS

df.pref.nap.ns = df.pref %>% 
  select(var_name, p_val_d_nap_ns) %>% 
  mutate(comparison = "Nap vs. NS") %>% 
  rename(p_val_unadjust = p_val_d_nap_ns) %>% 
  arrange(p_val_unadjust)

p.unadjust = df.pref.nap.ns$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.pref.nap.ns = df.pref.nap.ns %>% 
  bind_cols(p.adj.mat)

df.pref.all = df.pref.ns.cont %>% 
  bind_rows(df.pref.nap.cont) %>% 
  bind_rows(df.pref.nap.ns) %>% 
  relocate(var_name, comparison)

df.pref.all = df.pref.all %>% 
  rename(BKY06 = TSBH) %>% 
  arrange(comparison, var_name)

export(df.pref.all, "Datasets/pval_adj_pref.xlsx")

```


Correlation pvals

```{r}

#path = "C:/Users/pedro/Dropbox (Personal)/Sleep Data/01. Main Study/Analysis/final_paper/mht/Output"

df.pvals.sim = import("Datasets/mht_final_iv.dta")

df.nap = df.pvals.sim %>% 
  select(contains("_nap") )

cov.mat = cov(df.nap)
corr.mat = cov2cor(cov.mat)  

df.ns = df.pvals.sim %>% 
  select(contains("_ns") )

cov.mat = cov(df.ns)
corr.mat = cov2cor(cov.mat)  

indexes_nap = c("tstat_wellbeing_nap", "tstat_pref_nap", "tstat_cogfunction_nap", "tstat_earnings_nap")

df.ind.nap = df.nap %>% 
  select(all_of(indexes_nap))

cov.mat = cov(df.ind.nap)
corr.mat = cov2cor(cov.mat)  

indexes_ns = c("tstat_wellbeing_ns", "tstat_pref_ns", "tstat_cogfunction_ns", "tstat_earnings_ns")

df.ind.ns = df.ns %>% 
  select(all_of(indexes_ns))

cov.mat = cov(df.ind.ns)
corr.mat = cov2cor(cov.mat)  

# P values
df.pvals.only = df.pvals.sim %>% 
  select(p_wellbeing_d, p_pref_d, p_cogindex_d, p_earnings_d)

cov.mat = cov(df.pvals.only)
corr.mat = cov2cor(cov.mat)  


```

